/********************************************************************
* A new Lower Bound: Casiano
********************************************************************/
// Added benefit subject to constrains
unsigned int * cross(unsigned int * x, unsigned *y) {
  unsigned int * z = (unsigned int *) checked_malloc(sizeof(unsigned int)*pbm->n);

  unsigned int i;
  for(i=0; i<pbm->n; i++) {
    z[i] = x[i] + y[i];
    if (z[i] > pbm->x[i]) { 
      z[i] = pbm->x[i];
    }
  }
  return (z);
}

// Returns a fresh empty set
unsigned int * Init_Set() {
  unsigned int * set = (unsigned int *) checked_malloc(sizeof(unsigned int *) * pbm->n);
  unsigned int k;
  for(k=0; k < pbm->n; k++)
    set[k] = 0;
  return (set);
}

// returns the profit of a given solution
Bound profit(unsigned int * sol) {
  unsigned int i;
  Bound sum = 0;

  for(i = 0; i < pbm->n; i++) {
    sum += sol[i]*pbm->p[i];
  }
  return (sum);
}

#define NONE (pbm->n+1)

void CLB(void) {
  Bound **E; // Heuristic solution value for a metarectangle of dim x,y
  unsigned int ***S; // S[x][y] = Heuristic solution set for a metarectangle of dim x,y
  unsigned int x, y, k;
  
  E = (Bound **) checked_malloc(sizeof(Bound *) * (pbm->L + 1));
  S = (unsigned int ***) checked_malloc(sizeof(unsigned int **) * (pbm->L + 1));

  // Initialize E and S
  for (x = 0; x <= pbm->L; x++) {
    E[x] = (Bound *) checked_malloc(sizeof(Bound) * (pbm->W + 1));
    S[x] = (unsigned int **) checked_malloc(sizeof(unsigned int *) * (pbm->W + 1));
    Bound bestpiece;
    unsigned int k0;
    for (y = 0; y <= pbm->W; y++) {
      S[x][y] = Init_Set();
      bestpiece = 0;
      k0 = NONE;
      for (k = 0; k < pbm->n; k++)
        if ((pbm->l[k] <= x) && (pbm->w[k] <= y) && (pbm->p[k] > (unsigned int)(bestpiece))) {
          bestpiece = pbm->p[k];
          k0 = k;
        } /* end if */
      E[x][y] = bestpiece;
      if (k0 != NONE) S[x][y][k0] = 1;
    }
  }

  for (x = 0; x <= pbm->L; x++) {
    for (y = 0; y <= pbm->W; y++) {
      Bound aux;
      unsigned int * bs; // best solution
      Bound total;
      unsigned int max_x1 = x / 2;
      unsigned int max_y1 = y / 2;
      unsigned int x1;
      unsigned int y1;

      // Initial values
      aux = E[x][y];
      bs  = S[x][y];

      // Vertical
      unsigned int * combined;
      for (x1 = 1; x1 <= max_x1; x1++) {
        combined = cross(S[x1][y], S[x - x1][y]);
        total = profit(combined);
        if (total > aux) {
          aux = total;
          bs = combined;
        }
      } // end for

      // Horizontal
      for (y1 = 1; y1 <= max_y1; y1++) {
        combined = cross(S[x][y1], S[x][y - y1]);
        total = profit(combined);
        if (total > aux) {
          aux = total;
          bs = combined;
        }
      }

      // Best
      E[x][y] = aux;
      S[x][y] = bs;
    } // for y
  } // for x
  
  if (E[pbm->L][pbm->W] > bestSolution) { // Modify global variables 
    bestSolution = E[pbm->L][pbm->W];     // containing the initial solution
    printf("Heuristic = %ld\n",bestSolution);
  }

  // Free memory
  for (x = 0; x <= pbm->L; x++) {
    for (y = 0; y <= pbm->W; y++) {
      free(S[x][y]);
    }
    free(E[x]);
    free(S[x]);
  } // for x
  free(E);
  free(S);
}
